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We report on our progress in using the overlap-Dirac fermion operator in 
simulations of lattice QCD. We have investigated the Lanczos based method of 
Borici, as well as various rational approximations, to calculate the step function 
in the overlap-Dirac operator. The QCD simulations were performed at (5 = 
5.85 with a lattice volume of 8 3 32. 
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I. INTRODUCTION 

The results of lattice QCD calculations of weak matrix elements are a critical com- 
ponent of the experimental program in heavy flavour and kaon physics. The results from 
lattice gauge theory calculations would significantly improve if the masses of both the sea 
and valence quarks could be reduced. Unfortunately progress in doing this is very slow 
with simulations that use either the Wilson or clover fermion operators [l|]. 

It seems plausible that the difficulty of simulating with light quark masses with the 
clover and Wilson fermions operators is due to explicit chiral symmetry breaking in the 
actions. If the fermion operators were invariant under chiral symmetry transformations, 
their eigenvalue spectrum would be constrained to a smaller region || . The performance 
of the simulation algorithms degrades as the range of eigenvalues gets larger. Simulations 
that use the staggered fermion operator can reach much lower quark masses |I| than 
simulations that use the Wilson or clover operators, because the staggered action has a 
residual of the continuum chiral symmetry. Neuberger has derived || a fermion operator, 
called the overlap-Dirac operator, that has a lattice chiral symmetry [[|, [5[ . 

Our goal is to simulate QCD with the overlap-Dirac operator in the mass region: 
{Mps/My = 0.3 — 0.5). This quark mass region is inaccessible to lattice QCD simulations, 
that are currently computationally feasible, with the clover or Wilson operators |[|. 

* Talk presented at Chiral '99, Sept 13-18, 1999, Taipei, Taiwan. 
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II. THE OVERLAP-DIRAC OPERATOR 



The massive overlap-Dirac operator [|3], Q is 

^-id+P + d-^-T^—) (i) 

z \ H (m)H(m) 



where H(m) is the hermitian Wilson fermion operator with negative mass, defined by 

H(m) = -f 5 (D w - m) (2) 

where D w is the standard Wilson fermion operator. The parameter \x is related to the 
physical quark mass and lies in the range to 1. The m parameter is a regulating mass, 
in the range between a critical value and 2. 

III. NUMERICAL TECHNIQUES 

Quark propagators are calculated using a sparse matrix inversion algorithm. The 
inner step of the inverter is the application of the fermion matrix to a vector. For com- 
putations that use the overlap-Dirac, the step function 

<**-jfa (3) 

must be computed using some sparse matrix algorithm. A number of algorithms that 
calculate quark propagators from the overlap-Dirac operator, without using a nested al- 
gorithm have been proposed |7], ||. We discuss one of them in section [VI. 



Practical calculations of the overlap operator are necessarily approximate. To judge 
the accuracy of our approximate calculation we used the (GW) Ginsparg- Wilson error: 

1 1 l5 D N x + D N l5 x - 2D N j 5 D N x 1 1 j^-rr (4) 

|| x || 

which just checks that the matrix obeys the Ginsparg- Wilson relation H. Other groups P, 



Hj have found that the step function must be calculated very accurately, so we also use 
more sophisticated estimates of the numerical error (see Eq. |S[). 

Our numerical simulations were done using (3 = 5.85 quenched gauge configurations, 
with a volume of 8 3 32. This allows us to directly compare our results with the two other 



QCD spectroscopy calculations || jlOfl . The quark propagators were generated from point 



sources. For all the algorithms we investigated, we used m equal to 1.5. Although we are 
currently investigating using the overlap-Dirac operator in the quenched theory, most of 



the algorithms can also be used in full QCD simulations [11, [I 
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IV. LANCZOS BASED METHOD 



Borici |13| has developed a method to calculate the action of the overlap-Dirac 
operator on a vector, using the Lanczos algorithm. In exact arithmetic, the Lanczos 
algorithm generates an orthonormal set of vectors that tridiagonalises the matrix. 

HQ n = Q n T n (5) 

where T n is a tridiagonal matrix. The columns of Q n contain the Lanczos vectors. 

The "trick", to evaluate the step function (Eq. [3]), is to set the target vector b, as 
the first vector in the Lanczos sequence. An arbitrary function / of the matrix H acting 
on a vector is constructed using 

(f(H)b) t = ^Qnf{T n )Qi) l3 b 3 (6) 

= WblKQnfiTn))^ (7) 

where the orthogonality of the Lanczos vectors has been used. The f(T n ) matrix is com- 
puted using standard dense linear algebra routines. For the step function the eigenvalues 
of T n are replaced by their moduli. Eq. |7] is linear in the Lanczos vectors and thus can be 
computed in two passes. 

The major problem with the Lanczos procedure is the loss of the orthogonality of 
the sequence of vectors due to rounding errors. It is not clear how this lack of orthogo- 
nality effects the final results. Some theoretical analysis has been done on the use of the 
Lanczos algorithm to calculate functions of matrices Jl4|]. It is claimed that the lack of 



orthogonality is not important for some classes of functions. 

On small lattice we checked that the eigenvalue spectrum of the overlap-Dirac oper- 
ator moves closer to a circle [D|, as the number of Lanczos steps increases. Even after 50 



iterations of the Lanczos algorithm, there are still small deviations from the circle. For a 
hot gauge configuration with a volume of 4 4 , all the Lanczos vectors were stored and then 
were used to investigate the effect of the loss of orthogonality. The scalar product be- 
tween two Lanczos vectors drops from 10~ 7 to 10~ 3 after about 130 iterations, indicating 
problems with orthogonality. The GW error was 0.11 at 50 iterations and 0.012 at 250 
iterations. This is some evidence that the Borici's algorithm may still work even when the 
orthogonality of the Lanczos vectors is lost. It is much harder to look at the eigenvalue 
spectrum of the overlap-Dirac operator on a 8 3 32 (3 = 5.85 gauge configuration, so we 
computed the GW error instead. The GW error on a single gauge configuration was: 
2 10- 3 (100 iterations), 1 10~ 3 (200) iterations, 1 10" 4 (300) iterations, and 2 10~ 6 (500 
iterations). 
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FIG. 1. Effective mass plot for pion using Borici's method of calculating the step function using 
100 Lanczos iterations. 

Fig. |], is an effective mass plot of the pion, for four different values of the quark 
mass fi. The number of iterations in the Lanczos algorithm was kept constant at 100. 
The effective mass plots for the pion in Fig. [I] can be compared to the two other published 
spectroscopy calculations, that both use configurations with parameters: 8 3 16 and (3 = 



5.85. For a quark mass of \i = 0.1 (0.05), Liu et al. [|T(J obtain a pion mass in lattice units 
of 0.63 (0.45). From the graph by Edwards et al. ||, the pion mass was 0.87 (0.63) at 
/i = 0.1 (0.05). The differences between the two groups are probably explained by them 
using different values of m in their simulations, because the quark mass [i is renormalized 
by a multiplicative factor that depends on the domain mass ||. The effective mass plots 
in Fig. [I] are consistent with the data of Edwards et al. 0. Although smeared correlators 
should be used for a more detailed comparison. The quality of the \i = 0.03 effective mass 
plot of the pion is disappointing. The inversion of the overlap-Dirac operator that used 
100 Lanczos iterations, at a mass /i = 0.1, required 150 iterations in the inverter and took 
105 minutes on 32 nodes of our cray t3e. 

We checked the stability of the pion's effective mass with the number of iterations 
used in the Lanczos procedure. The pion effective mass plot was stable for the quark 
masses /i = 0.1 and 0.03, as the number of Lanczos iterations were varied from 100 to 
300. 



Liu et al. |H| used the Gell-Mann-Oakes-Renner (GOR) relation, that was derived 



in B for the overlap-Dirac operator, as a check on the accuracy of the computation of 
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FIG. 2. The ratio of the right hand side over the left hand side in equation || for a single 
configuration. The error bars are from the Z2 noise method used to compute the chiral 
condensate. The step function was computed using Borici's method with: 50, 100, 200 
iterations. 

the step function. 

MEW af M°)> = ^E^W(*)> ( 8 ) 

x V X 

where 7r is the pion interpolation field, and x is summed over the space-time volume (V). 
The "external" quark propagators |J denned by 

DM = j^P-'Oj) ~ 1] (9) 

should be used in equation 

The data in Fig. |2| show the GOR relation is satisfied up to 2% for the masses 
\i = 0.1 and 0.15, and 4% for the mass fi = 0.05. Increasing the number of Lanczos 
iterations did not decrease the violation of the GOR relation. This may be due to the 
loss of orthogonality in the Lanczos vectors. 
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V. RATIONAL APPROXIMATION 



The step function can be approximated by a rational approximation [JT6], I7 |. 

<tf)~tf(* + g^) do) 

The rational approximation is only an accurate approximation to the step function in a 
certain region. The eigenvalues of the matrix H should lie in this region. The coefficients 



Cfc and dk can be obtained from the Remez algorithm ||18|| . The number of iterations 
required in the inverter is controlled by the smallest dk coefficient, that acts like a mass. 
If dk is small, the number of iterations required is controlled by the condition number of 
H 2 . 

On one configuration we obtained GW errors of: 8 10~ 5 , 1 10~ 5 , and 2 10~ 6 , for 



the N = 6, N = 8, and N = 10, optimal rational approximations Unfortunately, 
the above results required up to 600 iterations for the smallest dk, which was too large to 
use as the inner step of a quark propagator inverter. We have not yet implemented the 



technique of projecting out some of the low lying eigenmodes ||18|| . This projection will 
reduce the condition number of the matrix, and hence the number of iterations required 
in the inner inversions. 

VI. Five dimensional representation of the overlap-Dirac operator 

One undesirable feature of the algorithms just presented for inverting the overlap- 
Dirac operator is that at each iteration of the inverter, some sparse matrix techniques 
must be done to calculate the overlap-Dirac operator. In the language of Krylov spaces, 
the theory of which underlies the numerical calculations, two independent Krylov spaces 
are used in a nested inverter. If the overlap-Dirac operator could be calculated using one 
Krylov space, this may be more efficient. Neuberger has proposed one method to calculate 
the overlap-Dirac operator without the nested inversion [|7|. The first implementation of 
Neubeger's ideas was discussed by Edwards at this meeting. 

To explain the idea we will use a simplified rational approximation. The generali- 
sation to higher order rational approximations is obvious. 

The equation for ip can be solved using additional variables (0j). The additional equations 
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generated by the new variables can be written in matrix form. 
/ 
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(12) 



where we have introduced the notation H c = H 2 + c. 

The additional variables makes the calculation of the overlap-Dirac operator look 
similar to the calculation of the domain wall operator [[19]. Although with an accurate 
enough rational approximation, this technique will calculate the overlap-Dirac operator 
exactly. The key issue is the condition number of the five dimensional matrix, because 
this controls the number of iterations required in the inverter. As the various rational 
approximations use small coefficients, these could have a large effect on the condition 
number. 

To study the effect of the rational approximation on the condition number of the 
five dimensional matrix, we have started to study the problem in free field theory. The 
calculation of the eigenvalues of the matrix in Eq. [12] is simple in free field, because Fourier 
analysis can be used. Although the free field theory eigenvalues will not be too similar to 
those of the interacting theory (although projecting out the lowest topological eigenmodes 
will improve the agreement), this is the only case where we have any hope of analytical 
insight into the condition number of the five dimensional matrix. 

The five dimensional matrix was constructed using MATLAB, using the free her- 
mitian Wilson operator [pD| . Table [I] contains some results for a 8 4 lattice, with a quark 
mass of /i = 0.1, and a domain mass of m = 1.0. The NIG approximation is the 16th order 
approximation to the step function introduced by Neuberger and Higham [01], [T7|. The 



R6, R8, and RIO rows are for the Remez approximations to the step function introduced 
by Edwards et al. | I8|| . The validity column is the maximum distance that the rational 
approximation deviates by 10~ 3 from unity, divided by the minimum distance. This is a 
measure of how good an approximation the rational function is. The order in Table [T] has 
been normalised so that it is comparable to the length of the lattice in the fifth dimension 
for domain wall fermions (the true order is obtained by multiplying by 12). 

Table [I] shows that the condition number of the five dimensional matrix strongly 
depends on the type of rational approximation used to construct it. It is interesting 
to compare the A^16 and R8 approximations that are almost equally good, but which 
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Approximation 


Validity 


Order 


Condition number 


N16 


69 


32 


790 


R6 


16 


13 


1460 


R8 


61 


17 


2580 


RIO 


303 


21 


4260 



TABLE I. Condition number of the five dimensional matrix 



have very different condition numbers. It would be instructive to compare the condition 
numbers in Table [l] with the condition number of the domain wall fermion operator [21 



VII. CONCLUSIONS 



The results for the masses of the light hadrons obtained by Liu et al. [[HJ and 



Edwards et al. || seem to show that QCD simulations can be run at much lighter quark 
masses, than can be explored with the standard Wilson or clover operators. To work 
with quarks as light as those in the calculations of Edwards et al. and Liu et al., we 
need to project out the lowest eigenvalues from the H matrix. We are working on an 
implementation of the eigenvalue projection technique. 

What is not clear is how expensive the simulations with the overlap-Dirac operator 
will be on more realistic lattice volumes. The only way the Wilson or clover operators can 
be used to simulate QCD with lighter quarks is by the brute force approach of simulating 
closer to the continuum limit. This too will be very expensive. As the overlap-Dirac 
operator has a lattice chiral symmetry, it should be able to be used to explore the light 
quark mass region of QCD in an elegant way. 

This work is supported by PPARC. The computations were carried out on the T3E 
at EPCC in Edinburgh. We thank K. Liu, T. Kennedy, R. Edwards, C. Michael, A. Irving 
and U. Heller for discussions. 
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